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We analyze the Hamiltonian time evolution of classical SU(2) Yang-Mills-Higgs theory with a 
fundamental Higgs doublet on a spacial lattice. In particular, we study energy transfer and equili- 
bration processes among the gauge and Higgs sectors, calculate the maximal Lyapunov exponents 
under randomized initial conditions in the weak-coupling regime, where one expects them to be 
related to the high-temperature plasmon damping rate, and investigate their energy and coupling 
dependence. We further examine finite-time and finite-size errors, study the impact of the Higgs 
fields on the instability of constant non-Abelian magnetic fields, and comment on the implications 
of our results for the thermalization properties of hot gauge fields in the presence of matter. 
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O ' I- INTRODUCTION 

(N : 

' A variety of essential physical processes, ranging from ultrarelativistic heavy-ion collisions fH to the reheating period 
. and phase transitions in the early Universe proceed at least initially far from thermodynamic equilibrium and 

involve abundantly many nonperturbative degrees of freedom. The first-principle based theoretical treatment of such 
, phenomena, which require a quantum field theoretic description but are inaccessible to Euclidean lattice simulations, 

is as a rule beyond present capabilities. Important exceptions to this rule arise, however, if the underlying amplitudes 
I— i] receive dominant contributions from classical fields. The latter may be provided, in particular, by bosonic long- 

• wavelength modes at high temperature T and with energies E <^ T since the Bose-Einstein distribution supplies 
' I them with the large occupation numbers needed to ensure (semi-) classical behavior. In non-Abelian gauge theories, 

observables governed by such classical modes are typically of O {g^T^ in the weak-coupling regime (where g is the 
' gauge coupling and g^T sets an inverse classical length scale) and have a finite classical limit. Prominent examples 
'— 'i include the transport coefficients which control magnetic screening Q and color diffusion [4], and in particular the 
static plasmon damping rate The latter has direct impact on the local energy and momentum equilibration 
processes among hot gauge-field quanta, which were found to occur over surprisingly short times of less than 1 fm/c 
in the excited matter created by ultrarelativistic nuclear collisions at RHIC P, 0] . 
' Essentially classical nonequilibrium observables of the above type may therefore be calculated by relating them to 
0^ real-time evolution properties of classical long-wavelength gauge fields and by simulating those nonperturbatively on 
' a spacial lattice Along these lines, the plasmon damping rate was argued to be proportional to the classical 

- gluon damping rate and, at least at weak coupling, further to the maximal Lyapunov exponent (MLE) which governs 
the exponential separation rate between initially neighboring random gauge-field configurations 0, Q . The underlying 
0^ , reasoning is based on the expected ergodicity of the classical field trajectories and on the relation between exponential 
' growth and damping rates provided by time reversal symmetry In the weak-coupling region, furthermore, these 
^ ' relations can be tested quantitatively by coinparison with results from partially resummed thermal perturbation 
• '~j theory or alternatively from kinetic theory . 

rN , Following up on the above arguments, the present paper will deal with the real-time evolution of classical SU(2) 

• Yang-Mills-Higgs (YMH) theory on spacial lattices of various sizes. A particular focus will be on the role of the scalar 
" " ' and hence classically treatable matter fields, provided by the fundamental Higgs doublet, in the chaotic dynamics. 

The center piece of the analysis is a systematic survey of the energy and coupling dependence of a set of maximal 
Lyapunov exponents designed to cover representative parts of the weakly coupled YMH phase space. Since our 
theory corresponds to the electroweak sector of the standard model with vanishing Weinberg angle, the resulting 
MLEs contain information which may be useful for understanding cosmological nonequilibrium processes during 
semiclassical evolution phases of the early Universe, including topological structure formation, baryogenesis [Tlj and 
potentially cosmic string evolution [T^ . 

Moreover, our results will be relevant for the analysis of local equilibration processes in the highly excited matter 
produced at the RHIC P, Q and soon the CERN LHC [is] colliders. Indeed, the chaoticity of the gauge dynamics 
provides a natural mechanism for entropy production by soft fields (and the accompanying particle production in the 
quantum case), and its most unstable field modes contribute dominantly to equilibration processes. In particular, 
our results will give rise to new estimates for the energy and coupling dependence of the gauge-field damping rate 
in the presence of scalar matter. Furthermore, the MLEs should receive contributions from the non-Abelian plasma 
instabilities which were recently argued to accelerate the isotropization and thermalization processes in the aftermath 
of high-energy nuclear collisions [14| . The underlying unstable modes could in principle be isolated by numerical 
techniques similar to ours. As in chaotic inflation scenarios, furthermore, such instabilities typically generate nonper- 
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turbatively large occupation numbers, which may extend the rehabihty of our classical treatment to larger couplings 
and lower temperatures. Some of our qualitative results may even be robust enough to provide guidance on the impact 
of fundamental quark fields. 

Although our main focus will be on the evolution of random fields, we also study the impact of the Higgs fields 
on the instability of a constant non-Abelian magnetic field. The employed techniques may later be applied to more 
complex coherent fields, including classical solutions of YMH theory [151] and gauge- invariant coherent soft modes p^ . 
Studies of this type could provide new insights into the corresponding quantum theories. Applied to multimonopole 
configurations of YMH theory with an adjoint Higgs field, whose chaotic interactions we have recently studied [17| . 
they may for example help to clarify the role of chaotic monopole ensembles in disordering the gauge-theory vacuum. 

The paper is organized as follows: in Sec. |TT]we summarize the formulation of SU(2) YMH theory on a Hamiltonian 
lattice, derive the corresponding field equations and discuss suitable distance measures on the space of gauge and 
Higgs field configurations. Section Hill outlines the main ingredients of our numerical analysis, examines finite-time and 
finite-size effects, discusses the time evolution of the energy transfer between the various field sectors, and evaluates 
the rate of divergence between initially neighboring random field configurations at intermediate times. On this basis, 
we generate in Sec. IIVI a representative set of maximal Lyapunov exponents, discuss their energy dependence and 
relation to the plasmon damping rate, then extend the analysis by calculating a set of long-time Lyapunov histories, 
and finally evaluate the impact of the Higgs fields on the Savvidy instability of constant non-Abelian magnetic fields. 
Sec. |V] puts our results into context by discussing related nonequlibrium processes in the early Universe and in the 
aftermath of high-energy nuclear collisions, and Sec. I VII summarizes our main findings and provides some conclusions. 



II. YANG-MILLS-HIGGS DYNAMICS ON A SPACIAL LATTICE 



In order to identify and measure chaotic properties of a dynamical system, one has to follow the evolution of its 
dynamical variables over sufficiently long periods of time. A numerical treatment of field theories further requires to 
approximate space by a discrete lattice. The analogous handling of the time variable (as typically implemented in 
Euclidean spacetime subject to periodic boundary conditions) is unsuitable for chaos investigations, however, since 
it would unacceptably restrict the accessible evolution times. Hence we resort to the Hamiltonian formulation of 
lattice field theory [l8| in Minkowski space where gauge fields are restricted to temporal gauge and time remains 
an unbounded and (in principle) continuous variable. A further benefit of this formulation is that residual gauge 
symmetries enforced by Gauss' law can be accurately preserved during time evolution. In the following subsections 
we briefly summarize this approach as it applies to YMH theory and define the distance measures needed to determine 
the Lyapunov exponents. (More details can be found e.g. in Refs. |J-&, ^1^.) 



A. Hamiltonian lattice setup 

In the following section we outline pertinent aspects of the Hamiltonian formulation of 3-f 1 dimensional SU(2) 
Yang-Mills-Higgs theory on a spacial cubic lattice subject to periodic boundary conditions. Since the Higgs field is 
taken to transform in the fundamental representation of the gauge group, this theory is equivalent to the electroweak 
sector of the standard model in the limit of vanishing Weinberg angle. The gauge is fixed to Aq = 0, i.e. to Weyl 
gauge. The unbroken phase corresponding to a gauge-matter plasma is selected by positive Higgs mass and interaction 
terms, which allows for comparison of the results with (hard-thermal-loop resummed) perturbative results at high 
temperature below. 

The corresponding YMH Hamiltonian can thus be written as 

x;i X l<i<j<3 

X x,i 

+a ^ ^tr (0t 0^) [6 + a^^ltr {4>l4>x)] (1) 

X 

where g is the gauge coupling, k, the Higgs self-coupling, a the lattice spacing and dots denote time derivatives. The 
non-Abelian magnetic field is described by the spacial plaquette 

U^,^3 = U,,^ U^+^,, U-l^^^ f/-j ^ exp {-iga^F^^^j + O (a^) ) (2) 
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(with i,j € {1, 2, 3} and i ^ j), i.e. by the ordered minimal-circumference loop constructed from the link variables 

Ux,i = cxp {-igaA^.i) (3) 

where A^^i — A"^ ^t°- is the gauge field and cr" = with a G {1, 2, 3} are the Pauli matrices. The Ux,i are defined on 
the link which connects the site x with its neighbor in the positive i direction. Hence the spacial plaquettes contain 
the non-Abelian magnetic field strength components Fx^ij = \ {A^+ij — A^j ~ A^+j.i + A^^i) — ig [A^^i, A^j] while 
their electric counterparts ^ — —A°^ ^ are independent variables. The first term of the Hamiltonian ([T|) therefore 
describes the energy residing in the electric fields while the second term, 

X l<i<j<3 x,i,j 

approaches the magnetic or potential energy of the gauge field in the naive continuum limit. 

For the numerical implementation of the SU(2) link variables we have adopted the quaternion representation 

jj - ( u° -iu^, -u^ -iu^\ , , 

(7 = w — • CT = 2 • 1 , • 3 5) 
y — lu'-, vr + lu'^ J ^ ' 

(the indices x,i are suppressed) whose real components — (u'^,uj ^ R, fi £ {0,1,2,3} satisfy the constraint 
det U = M^u" + u"u" = 1 and thereby ensure unitarity UW = 1 as well. The u'^ are thus (four dimensional, 
cartesian) coordinates on the SU(2) group manifold S^. The representation ([5]) leads to simple field equations (cf. 
Sec. Ill B[) and requires the minimal number of floating point operations to calculate the product UV — — u°'v'^ — 
ia'^ (vPv'^ + v^u"' + e°'^''u^v'^) of two link variables. In order to state the initial conditions for the time evolution of 
the gauge field, however, we prefer the alternative representation of the link variable as a rotation of angle ujq around 
the direction n (d, ip), i.e. 

, ( — IT) ■ /T Sin ( 1 



U = exp [^igA^—j = cos (^— j - in ■ asin j (6) 

(suppressing again the indices). In terms of the polar angles < loq < 2tt, < i9g < ^ and < ipG < 2tt one then 
has gA'^ = uiqu'^ with n°' = (siniSc cos ipQ, sini^Q sin ipQ, cosi?g) and m° = cos {ujq/2), = sin {ujq/2). The Higgs 
field in the fundamental representation of the gauge group is written in an analogous quaternion representation, 



R 



cos — zn • cr sm 

\ 2 J V 2 



(7) 



where the polar decomposition again turns out to be more suitable for stating the initial conditions (cf. Sec. IIII A[) . 
In contrast to the unitary link variables U, however, the (square) modulus 

i?2 = ^tr (0t0) (8) 

of the Higgs field remains unconstrained. 

Exploiting its (classical) scaling properties, the YMH Hamiltonian ([1]) can be reexpressed in terms of the dimen- 
sionless variables H = g^aH, E% ^ — ga^E^^, (f>x — ga(f>x, k = n/g'^ and i = t/a as 

-ff = E [^G.cl (x) + £G,inag (x) + eH,kin (x) + EH.pot (x) + EG-H (x)] (9) 

X 

with the dimcnsionless energy densities 

ecoiix) ^ eG,mag(a:)-4 ^ fi-ltrUx,^X (10) 

i l<2<i<3^ ^ 

eH,kin(a;) = ^tr ($l$x^ , EG-nix) = -^tr {(j>l.Ux.i4>x+i) , (11) 

i 

EH.pot (x) = 3tr {4>l4>x) + J [tr {^Ux)] ' , (12) 

where the fields are now functions of t and dots represent d/ dt. The above form of the Hamiltonian renders the 
dependence on the total energy H and the Higgs self-coupling R, i.e. the two physical parameters of the YMH system, 
explicit (whereas the lattice spacing a and the gauge coupling g are absorbed into the dimcnsionless variables and 
fields). 
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B. Field equations 

The YMH Hamiltonian ([9]) generates the classical time evolution of electric, magnetic and Higgs fields. This becomes 
explicit in the corresponding first-order Hamilton equations which we derive with the help of the Poisson brackets 

(13) 

dqs dps dps dqs 

of the dynamical variables X with the Hamiltonian H (where qsTPs are the canonically conjugate variables and 
summation over s is implied). According to the canonical formalism, the time dependence of X is then determined 
by its Hamilton equation 

X^\{X,H]. (14) 
9 

Specializing Eq. p4p to the link variable Ux,i and abbreviating I = {x, i} leads with 

{Ef,U^} ^-igH'^U^Si^ (15) 

to the equation of motion 

iyi^\{Ui,H} =iEiUi (16) 
9 

where Ei — Eft"". In the quaternion representation ([5]) this equation reads 

^? = \Etut, < = - i (Efuf + e'^'^E^ut) (17) 
and maintains, in particular, the time-independence of the miitarity constraint, i.e. 

u'}u^i+u^iU^ = Q. (18) 

Hamilton's equation for the non-Abelian electric field strengths E% j, which are the canonically conjugate momenta 
of the link variables, similarly becomes 

k., = \ {K.,^H] = \Y.^r [a^ - U-l)] + Ut {^1<7^U,,,^,+,) (19) 

where the sum goes over the four plaquettes which contain the link {x, i}. 

The Hamilton equations for the Higgs field, its canonical momentum tt and their Hcrmitian conjugates are analo- 
gously found to be 



4 = ^ {^.,H} = ^Trt = iir (4) (20) 
$1 = ^{^lH} = in, (21) 



and 



'^l^x)] - 2 ^ Ux,i(t>x+i I 



(22) 



In order to prepare for an efficient numerical solution of this system, we rewrite it in terms of two second-order 
equations. 
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0^ = ~Ktr ((f>l(f>x) (t)^ - 6(t)x + 2^Ux,i(f'x+i, (24) 

i 

$1 = ~fitr{^l^x)^l-6 4'l + 2J24>l-^Ux-^,^. (25) 

i 

and then combine those, by addmg the Hcrniitian conjugate of Eq. ([25|l to Eq. (|24|) . into 

0^ = -Ktr (0^.(^0:) - 6 0a; + ^ ( + C^~ij j02;-i) • (26) 

i 

Finally, we recall that the full YMH dynamics in temporal gauge is only recovered after supplementing Hamilton's 
equations ([TT]), and (^5]) by Gauss' law 



3 r 1 



1=1 



(27) 



which acts as a constraint. Since its Poisson bracket with the Hamiltonian ([9]) vanishes. Gauss' law is preserved under 
time evolution. Above, we have defined the dimensionless non-Abelian charge density 



Im 



[kha^-'')] (28) 



carried by the Higgs field. 



C. Distance measures for gauge and Higgs field configurations 

The chaotic behavior of dynamical systems reveals itself in an exponential sensitivity of their time evolution to small 
changes in the initial conditions. The quantitative characterization of this sensitivity requires a distance measure on 
the field configuration space (i.e. a metric). More specifically, in the YMH system one has to monitor the separation 
between a reference configuration [Ui, cf)) and its neighbor ([//, cf)') — [Ui +SU,(j) + 6(j)). We will use individual distance 
measures in the gauge and Higgs sectors for this purpose, in order to determine the distance growth rate between two 
initially nearby gauge and Higgs field configurations individually. 

In the gauge sector, we adopt the gauge-invariant metric 

[Ui = ^ E I ^^^P - ^'^K I (29) 
P p 

(where A^p — is the total number of plaquettes on a lattice with A^ sites per spacial dimension) to measure the 
distance between gauge- field configurations. In the continuum limit the distance measured by the metric (|29p becomes 
proportional to the difference between the potential energies of reference and neighboring gauge fields. In the Higgs 
sector we employ the metric 

ciH[0,0']-^5] |(i?.)'-(i?;)'| (30) 

X 

which is gauge invariant as well. 

Since the lattice gauge group, and consequently the 3 (A"^ — l) A^'^ dimensional space of magnetic SU(A'c) gauge- 
field configurations on a lattice with A^ sites per dimension, is compact and of nontrivial topology, more and more 
field configurations approach the same distance dc when A^ increases. For the same reason, the distance ((29| is 
bounded from above, and for fixed total energy an analogous bound applies to the whole phase space. These bounds 
lead to an eventual saturation of the distance growth. Although this does not limit the principal effectiveness of the 
measures for determining the Lyapunov exponents (see below), it adds to the typical "finite-time" uncertainties 
encountered in their numerical analysis. Other sources of finite-time errors arise from the need to extrapolate the 
numerical results to the t ^ oo limit in which the MLEs are formally defined, and for N ^ oo from the exponentially 
growing distances between chaotic trajectories which eventually overburden the fioating point number representation 
capacities of any computer. 
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The standard approach for keeping finite-time errors of MLEs under control is to periodically rescale the distances 
[2l| after time intervals r and to extrapolate the numerical results for In dg h it) /t to infinite times. This approach 
has been used to calculate several MLEs in non-Abelian gauge theories [20. [23. 125L [2^ and to determine the whole 
Lyapunov spectrum on small lattices ^2^1 . We have adopted the same technique for the calculation of several long-time 
trajectories to be discussed in Sees. IIII Bl and IIV CI In these cases, we found it advantageous to employ the alternative 
distances measure 

d';^''\Ui,Er,Ul,E',] = IJ2 [{Ei^E'y + {Ui^Ulf][ (31) 

in the phase space of the gauge fields, which is a variant of the measure used in Ref. ^IQ], and 

1/2 



EE ( 



2 



(32) 



adopted from Ref. [26], in the Higgs sector. Of course, the resulting Lyapunov exponents should not depend on the 
choice of distance measure. We have checked this for several examples and confirmed that the deviations between the 
MLE values obtained from the metrics (|29|). ([30|) and ([31]) . (l32|) indeed remain well below the one-percent level. 

Nevertheless, even under rescaling the practically achievable evolution times remain limited by the available com- 
puter resources. In fact, even in pure SU(2) Yang-Mills theory (25j systematic extrapolation errors turned out to 
become negligible only after evolution times of the order of 10^ lattice units. To make matters worse, we will find 
below that the equilibration between gauge and Higgs fields proceeds at a far slower pace than among the gauge fields 
alone (cf. Sec. IIII Bp . and that as a consequence substantially longer evolution times are required to suppress such 
extrapolation errors in YMH theory. Adherence to one of our main goals, namely to calculate a rather exhaustive set 
of MLEs in the weak-coupling parameter and phase space, will therefore require a compromise. Indeed, to cover the 
relevant initial parameter space (on lattices of several different sizes) requires the calculation of O (lO^) trajectory 
pairs and thus forces us to limit the individual evolution times. 

Fortunately, size and systematics of finite-time errors can be estimated on the basis of the long-time energy balance 
(cf. Sec. IIII Bp and a few long-time orbits (cf. Sec. IIV Cp . Since we are mainly interested in the systematic energy- 
, coupling- and lattice-size dependence of the MLEs (rather than in their precise numerical values), furthermore, 
the competing goals of error suppression and calculability can be reconciled reasonably well. Our compromise will 
be to follow the majority of our distance histories only until they have saturated (without rescaling), which yields 
sufficiently accurate MLE estimates for most of our purposes. Rescaling will be used, on the other hand, for the 
long-time trajectories which we need to examine the energy transfer and equilibration processes between the gauge 
and Higgs field sectors in Sec. IIII B I and for the analysis of the MLE's finite-time errors and saturation properties in 
Sec. UFO 



III. FIELD INITIALIZATION, ENERGY BALANCE AND DISTANCE EVOLUTION 

In the following section we discuss in turn the initialization of the neighboring field configurations, the distribution 
of the total energy over the different field sectors, and the time evolution of the distance between initially adjacent 
field configurations. 



A. Initial conditions 



Our first task will be to generate a representative set of phase-space trajectories for pairs of specifically initialized 
reference YMH fields (C//,0) and their neighboring configurations {U[,4>') = [Ui + 5U,(f> + 54)) at sufficiently small 
distances dG[Ui,U'i\ and c?h[0, 0']- The resulting distance evolution histories will provide one of the foundations 
for our subsequent analysis of the maximal Lyapunov exponents. In the present subsection, we select a set of 77 
initial conditions for the reference trajectories such that the weak-coupling region of the YMH phase space is covered 
with sufficient resolution. In order to allow for direct comparison with a previously calculated MLE, we follow the 
initialization procedure of Ref. [26| . The resulting sample of field-pair trajectories will be considerably larger than 
that of preceding MLE calculations in gauge theories and include results from substantially larger lattice volumes 
(with up to 30^ sites). 
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In order to satisfy Gauss' law (|27|) initially (and consequently over the whole time evolution), we set the non-Abelian 
electric field and the time derivative of the Higgs field at the initial time — equal to zero, i.e. 

E-,,, (0) = 0, 4 (0) = 0, (33) 

which implies (0) = (cf. Eq. ([^5]) '). Hence the initial kinetic energies of all fields vanish while the potential energies 
are finite and ensure that the system starts far from equilibrium. The link variables Ui are initialized by randomly 
choosing the isospin directions tJgj'/'g of the gauge potential Af = LUG,if^fi'^G,h fG,i) from their full domains, while 
the initial value of the amplitude luq is chosen randomly over the restricted domain wq,/ S [0, 27r(5] with 6 < 1. The 
value of the parameter 6 therefore controls the average gauge-field energy per plaquette, Ep {S), which grows as 6^ 
for (5^1 and saturates in the limit (5 ^ 1 at the value Ep — 4 (cf. Fig. [1]). The upper bound on Ep arises 
from the fact that the magnetic part (j4|) of the Hamiltonian ([1]) is uniformly bounded by the SU(2) group volume, 
-f^mag = ag^H^ag < 24:N^. The Higgs field 0, finally, is initialized by choosing its angular variables ujh, z?h and ipH 
randomly from their full domains while keeping the dimensionless amplitude Rx = gaRx fixed at the same value R 
for all X. As a consequence, the initial (potential) energy of the Higgs field is determined by the amplitude R and the 
coupling K. 

The above initialization scheme characterizes any phase-space trajectory on a given lattice by three parameters 5, 
R and R which determine the average initial energy of both the gauge and Higgs fields. In addition, we will vary 
the lattice size, specified by the number N of sites per dimension, so that each of our field-pair histories can be 
uniquely labeled by a quadruple of values for 5, R, R and N . Our maximal lattice size with TV = 30 is chosen to 
substantially reduce potential finite-size effects of previous studies [13, [H, [l^l on considerably smaller lattices. The 
main benefit of the random angle initialization is that it equips the initial configurations with a specific average energy 
density, or equivalently with a temperature T which the fields will reach after equilibration. In our context, this is 
important because the temperature dependence of the MLEs is used to relate them to the static plasmon damping 
rates. Moreover, the resulting MLE values will turn out to be (within errors) independent of the random part of a 
given starting configuration, which indicates that the autocorrelation functions of the fields have decayed sufficiently 
strongly before the MLEs are measured (see below). 

In order to stay safely inside the validity range of the semiclassical approximation, and to be able to relate our 
findings to perturbative results, we will restrict our simulations to the weak-coupling regime. As pointed out in Ref. 
[1^ , this requires that the energy contributed by the Higgs mass term dominates over the Higgs self-interaction energy, 
i.e. 

< 1, (34) 

and that the magnetic gauge-field energy dominates over the gauge-Higgs interaction energy (which implies a weak 
gauge- Higgs couphng), or 

R^ <5 (35) 

(for maximal field amplitudes). Both conditions also improve the eventual equipartition of the electric, magnetic and 
Higgs field energies because they prevent the total energy from strongly exceeding the bounded magnetic energy. The 
lower bound (|35p on 5 additionally limits finite-size effects [cf. Eq. ([M)) ]. Furthermore, 5 values too close to 1 should 
be avoided in order to keep lattice-spacing artefacts under control and to remain sufficiently close to the continuum 
limit (cf. Eq. dSHD). 

For each of the initial reference configurations {Ui, (p) created according to the above procedure, we also generate a 
neighboring configuration {Ul,(j)') separated from {Ui,<f)) by distances dG{t = 0) < 5 x 10^'' and dB_{t = 0) < 10^^''. 
This is achieved by randomly choosing slight variations of all the reference configuration's field angles in the range 
5ujG, SSq, S(fG, Sujh, SSh, Sifiii £ [— e,e] where e — 10~^. We then integrate the field equations of Sec. Ill Bl for each 
of these configuration pairs [81| by means of a fourth-order Runge-Kutta algorithm and determine the time evolution 
of the distances c?g and dn- The integration time step should be much shorter than the lattice spacing a, and is 
additionally chosen small enough to ensure energy conservation with an accuracy of more than eight significant digits 
(after each step). The maximal violation of the constraints detJ7;= 1 after a single time step of length At = 10~^a 
is then about 10~^^ at each link. In order to avoid the accumulation of these round-off errors, we further rescale 
the link variables after each step such that their determinant remains exactly unity (and Eq. psp exactly satisfied). 
We convinced ourselves that Gauss' law (|27p then remains satisfied to better than five significant digits after each 
integration step. 
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B. Energy distribution over gauge and Higgs fields 



^^G — Eel + E^aag — -Ep, (36) 



We now turn to the energy transfer processes between the electric, magnetic and Higgs fields which contain crucial 
information on the nonequilibrium dynamics and quantitative thermalization properties of the YMH system. In our 
context, this information will be particularly helpful for understanding, qualitatively estimating and reducing the 
finite-time errors which afflict the calculation of the MLEs, and for putting the relation between the MLEs and the 
plasmon damping rates on a more solid footing. For several long-time trajectories, we have therefore recorded the 
evolution of the energies per degree of freedom stored in the electric field, E^i = Hd/ (6N^), in the magnetic gauge 
field, E'niag — -ffmag/ (SiV'^) , and in the Higgs field. En = Hh/ (4A^'^) , over the unprecedentedly long time periods 
< t < 20000a . In the following, wc will often express the total energy (per degree of freedom) Eq of the gauge field 
in terms of the average energy per plaquette Ep as 

1 

r 

and frequently encounter the total YMH energy per degree of freedom, E = {6Eq + 4£'h) /lO, as well. 

Typical results for the time evolution of the different energies are plotted in Fig. [5] with 5 = 0.2, R = 0.2, k — 1, 
N = 10 and in Fig. [3] with 6 = 1 and otherwise unchanged initial values. They confirm and extend the observation of 
Ref. [1^ that the energy equilibration between the electric, magnetic and Higgs field sectors of YMH theory proceed 
over two drastically different time scales (at least in the weak-coupling regime). Indeed, even when initialized in 
highly nonequilibrium configurations, as selected in Sec. IIII Al the electric and magnetic gauge sectors can be seen to 
equilibrate very rapidly, namely after only a few lattice time units a. The Higgs field's potential and kinetic energies, 
which are not shown separately in Figs. [2]and[3l equilibrate over an approximately equal relaxation time. (Generally 
the gauge and Higgs sectors reach different temperatures, however, according to the amount of energy stored in them 
by the initial conditions.) In contrast, the mutual thermalization of gauge and Higgs sectors typically requires 4 to 5 
orders of magnitude more time. In fact, the energy transfer between the two sectors becomes appreciable only after 
a few hundred time units and takes several thousand more to essentially complete for 6 = 1, and many more for 
6 = 0.2. Moreover, for the maximal i5 = 1 moderate deviations from complete equipartition of the energy remain 
visible in Fig. [3] even after 10000 time units have elapsed. This may be a consequence of lattice-spacing artefacts 
which are maximal a.t 6 = 1 (cf. Sec. IIII Ap . The huge discrepancy between the two characteristic relaxation scales 
can be largely attributed to the initial conditions of Sec. IIII Al which keep the system close to the weak-coupling and 
continuum limits. 

The gauge-field energy ([55]) can be directly related to the temperature T which the gauge fields reach after times 
t 3> A""'^ (where A is the MLE). At sufficiently weak coupling (among the field oscillators) one has j28|| 

for the gauge group SU(A^c) and thus Eq = T for = 2. This relation will be relevant for the evaluation and 
interpretation of the MLEs which we extract in Sec. IIV Al from the distance growth rates after the gauge fields 
became members of a prethermal ensemble. As mentioned in the introduction, lc\ = {g^T^ = (s^^'g) acts as a 
classical length scale in hot quantum gauge-theory amplitudes which depend (to leading order in thermal perturbation 
theory) on g and T exclusively in the combination o^T. This observation suggests additional conditions for keeping 
lattice artefacts in such amplitudes under control [2{|. More specifically, in order to remain sufficiently close to the 
continuum limit the lattice spacing a should be much smaller than ^ci, i-e. 

= ag^Ec « 1, (38) 

and in order to avoid finite-size effects the extent Na of the cubic lattice has to be much larger than Id, i.e. 

iV^G > 1- (39) 

As expected, these conditions require 3> 1, and the upper bound (j38p on Eq furthermore ensures that the underlying 
lattice structure cannot be resolved by the gauge fields. (Of course, for a ^ one will eventually encounter UV 
singularities of Rayleigh- Jeans type in some amplitudes, signalling the onset of indispensable quantum corrections to 
the classical field statistics.) 



C. Divergence of neighboring field trajectories in phase space 



In the following section we analyze the time evolution of the distances (Iq and djj between pairs of initially adjacent 
random field configurations which were generated according to the procedure of Sec. IIII Al and followed until saturation. 
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The MLEs and their parameter and in particular energy dependence will then be extracted from the growth rate of 
the logarithmic distances in Sec. IIV Al In order to cover the relevant phase space, we select a representative set 
of values for the parameters S, R, R and N which characterize any initial configuration. The initial, homogeneous 
Higgs field amplitude is fixed at i? = 0.2 for all trajectory pairs, which allows for a quantitative comparison with a 
configuration studied in Ref. [26| . To stay sufficiently close to the weak-coupling and continuum limits then requires, 
according to Eq. (f34|) . that the Higgs self-coupling is bounded by k < 25, and as a consequence of Eq. (f35|) that 
the initial magnetic (and total) gauge-field energy is restricted hy 6 > 0.04. As mentioned above, the bound on 5 
also helps to avoid significant finite-size artefacts [cf. Eq. and allows to extract the approximate MLEs with 

reasonable accuracy even after rather small evolution times (see below). 

Guided by the above arguments, we generate trajectory pairs for 11 values of (5 S [0.05, 1]. For each of them, we 
plot the resulting In c?g {t) (black lines) and In dn {t) (grey lines) in Fig. |3] at fixed Higgs self-coupling k = 1 on 
lattices of four different sizes corresponding to iV = 6, 10, 20 and 30, and in Fig. [5] on a TV = 20 lattice with the Higgs 
coupling values R — 1,8, 16 and 24. (The configuration pair studied in Ref. |,26j on a relatively small lattice with 
TV = 10 and 6 — 0.2, R — 0.2, k = 1 is therefore included in our sample.) The corresponding logarithmic distances 
IndcH {t) for the 11 5 values are grouped into three sets which are separately plotted in panels (a) - (c) of Figs. [4] 
andO in panel (a) we display lnc?G,H for the five largest values 5 = 1.0, 0.5, 0.45, 0.4, 0.35, in panel (b) for the values 
5 — 0.3,0.25,0.2,0.15 and in panel (c) for the two smallest values 5 = 0.1,0.05. All 5 values except for the smallest 
(i.e. 5 = 0.05, which is most strongly affected by finite-size artefacts) store more energy in the gauge than in the 
Higgs sector. 

The essential characteristic which all logarithmic distance histories of Figs. S] and E] share is that after a latency 
period of varying length they start to rise at least approximately linearly with t/a before reaching a time- independent 
saturation plateau (which lies somewhat outside the plotted i domain for 5 = 0.05) at the maximal distance in the 
compact phase space. Distance saturation at large t = t/a is a, consequence of the compactness of the lattice gauge 
group and could be avoided by periodical rescaling (cf. Sec. Ill C|) . The linear regions and the underlying exponential 
growth rates between initially almost identical field configurations reveal an exponential sensitivity of the distance 
evolution to the initial conditions, i.e. the standard hallmark of temporal chaos. Not surprisingly, the fields grow apart 
at a faster pace when their energy increases, i.e. the slopes in Figs. H] and [5] grow with &. For each field trajectory, 
furthermore, the linear regions of both In cLq (t) and In dfj (t) have the same average slopes. This result differs from a 
previous estimate for one trajectory 27] and will be discussed further in Sec. IIV Al Moreover, for 5 < 0.2 the latency 
period, which is hardly noticeable for larger 6, expands and the linear growth becomes increasingly modulated by 
oscillations whose frequency increases with S. This behavior was observed in YM theory as well and can be traced 
to the impact of the next-to-maximal Lyapunov exponents which grows when the maximal exponent decreases . 
Obviously, these oscillations reduce the accuracy with which the maximal Lyapunov exponent can be determined from 
the slopes of Inrfc.H (t) in the linear regions (see below). 

Figures m and [5] further show that for all field trajectories (except that with 6 = 0.05) Indn (t) stays below Indc (0- 
This reflects the smaller amount of energy initially stored in the Higgs sector for S > 0.05 (cf. Sec. IIII A[) and will 
change during the long-time evolution to be discussed in Sec. IIV CI In addition, the height of the saturation plateaus 
of lnc?G (t) decreases slightly with S while that of Indu (t) remains constant. This may indicate that the maximal 
magnetic gauge-field distance is reached only when sufficient gauge-field energy is available. Apart from these 
differences in the saturation behavior, however, even the modulation patterns of In (Iq {t) and In dfj (t) are very similar. 
This suggests that, despite the small gauge-Higgs coupling ensured by Eq. (|55|) . the time dependence of the gauge 
and Higgs components of at least the most unstable mode has already synchronized after a few lattice time units. 

The qualitative dependence of the results on the lattice size, i.e. on G {6, 10, 20, 30} with the lattice UV cutoff 
kept fixed, can be judged by comparing the distance histories in Fig. 31 Figure 3] a contains the results for 
I > S > 0.35. Although the fields are randomly initialized, the curves with identical 6 but different N clearly cluster, 
i.e. in accord with the bound ([55)1 essentially no finite-size effects can be observed in the covered N and 6 regions 
(while lattice-spacing effects should become noticeable for 6 close to unity [cf. Eq. (|38p]). Indications for a similar 
N independence were found in pure YM theory f^. This may suggest that the most unstable modes, i.e. those 
which dominantly drive the chaotic time evolution of initially adjacent configurations, have for sufficiently large initial 
magnetic field energy (corresponding to 6 > 0.35) typical wavelengths which are small enough to be accommodated 
by even the largest considered IR cutoff (corresponding to iV = 6), or in other words that these most chaotic modes 
essentially fit inside a periodic (6a) lattice volume. As shown in Figs. [4)d and |4|:, however, for smaller S < 0.2 
finite-size corrections become visible in the average slopes of the (increasingly oscillation modulated) linear regions of 
IndcH {t) and cause them to differ more strongly. A systematic trend in the N dependence of these slopes cannot be 



discerned in our data, however, whereas in pure YM theory the slope was found to increase on smaller lattices [25 1. 

Figure [5] reveals the qualitative dependence of the distance histories on the Higgs self-coupling R. In the range 
R G {1,8, 16,24} (for TV = 20) it bears several qualitative similarities with the TV dependence of Fig. [H To begin 
with, a R dependence is hardly noticeable for large S while the slope of the linear regions becomes increasingly R 
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dependent towards smaller values of 5, although again without a perceivable systematic trend. The, as a whole, only 
mild sensitivity of the slopes to R is probably a consequence of the fact that even In dn (t) is mainly determined by the 
most unstable gauge-field fluctuations and hence relatively insensitive to the self-interactions of the Higgs field. After 
full equilibration between the gauge and Higgs sector has taken place, the R dependence of the slopes may therefore 
be systematically enhanced (if distance saturation is avoided by periodical rescaling), as we will indeed find in Sec. 
IIV CI Since a more strongly self-coupled Higgs sector would absorb energy from the gauge sector (in which for S > 0.1 
more initial energy is stored, cf. Sec. IIII A[) faster, it should similarly increase the R dependence of the slopes. 

To summarize, all members of the representative set of distance histories (in the weakly coupled, symmetric YMH 
phase) discussed above increase exponentially and thereby exhibit chaotic behavior. For S > 0.3 the Indc.H (t) become 
practically independent of the Higgs coupling R and (for > 6) of the lattice volume. 



In the following section we proceed to the quantitative evaluation of the maximal Lyapunov exponents for random- 
ized and coherent initial conditions, and we discuss their energy dependence and relation to the plasmon damping 
rates. 



The analysis of the last section showed that all our 77 randomly initialized field pairs belong to the chaotic part 
of the YMH phase space. This suggests that in the unbroken phase of YMH theory chaotic behavior is either 
universal (i.e. exists for all energies) or at least prevalent in most of the weakly coupled phase space [s^. In order 
to quantify this behavior, we will now evaluate the classic measure for the chaoticity of a dynamical system, i.e. 
the maximal Lyapunov exponent A or equivalently the exponential growth rate of the distance between initially 
neighboring dynamical variables. We are going to extract the MLEs from the numerical results of Sec. Ill CI by 
averaging the time histories do {E, R, N; t) and c?h {E, k, N; t) of the gauge and Higgs field distance measures over 
the time interval A during which they remain in the linear regime, i.e. 



Note that we have replaced the dependence on the initialization parameter 5 with that on the total (dimensionless) 
energy E of the YMH system, and that we suppressed the dependence on the remaining initialization parameter, the 
Higgs amplitude R, which is kept at the same value for all our trajectories (cf. Sec. IIII A|) . We further recall that the 
above method for obtaining the MLEs becomes increasingly error-prone towards lower energies where equilibration 
proceeds more slowly while the impact of the next-to-maximal Lyapunov exponents grows and generates modulations 
of IndcH (t) with decreasing frequency. Similar problems were encountered in Ref. 7] and will be tamed below by 
periodically rescaling the distance measures (cf. Sec. IIV CI) . 

Our numerical results for the dimensionless MLEs Aq.h {E,R,N^ := aAcH {E,R,N^, based on the 77 field-pair 
evolution histories of Sec. IIII CI are collected in Table 1. A first glance at the table confirms the qualitative trends 
which we noticed in our discussion of Figs. |3]and[5]in Sec. IIII CI Besides the expected increase of the Ag,h with E (or 
(5), which we will analyze quantitatively in Sec. lIVBj the data show fiuctuations in the statistically expected range 
of about 10% for different Higgs self-couplings and lattice sizes, but except for the smallest E no obvious systematic 
dependence on either R or A^. At the considered intermediate times (i.e. after separate preequilibration of gauge 
and Higgs sectors but before their mutual thcrmalization is complete) and at least at intermediate energies E or 
Ep systematic finite-size and lattice-spacing effects are therefore small. Furthermore, the above results indicate that 
the Higgs sector plays a rather minor role in the chaoticity of the full YMH system, at least at the weak couplings 
which the initial conditions of Sec. IIII Al implement. The most unstable mode, which in large part drives the chaotic 
behavior, should therefore be controlled mainly by the gauge dynamics. As a consequence, reasonable estimates for 
the MLEs can be extracted at the preequilibration stage and the MLE values of the SU(2) YMH system should be 
similar to those of pure SU(2) YM theory 0,111, which is confirmed by the results in Table 1. 

While we find the Higgs sector to have only limited impact on the chaotic YMH dynamics in the weakly coupled 
symmetric phase, it may be useful to recall the results of Refs. [s^lssl] in the homogeneous limit, i.e. for wavelengths 
much larger than the inverse amplitudes \A\~ ,\4>\~ , which reveal a more dramatic role of the Higgs field in the broken 
phase (even at nonzero Weinberg angle). This is a consequence of the dynamically generated gauge-field mass in the 
broken phase which is known to damp (and beyond a critical value to fully suppress) chaotic behavior [20| . (We note in 



IV. LYAPUNOV EXPONENTS, SCALING BEHAVIOR AND DAMPING RATES 



A. Maximal Lyapunov exponents of randomly initialized fields 




(40) 
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iV = 6 
K = 1 


TV = 10 
K = 1 


N ^20 
K = 1 


TV = 30 
K = 1 


N ^20 
K = 8 


N ^20 
K = 16 


N ^20 
K = 24 


E = 0.04374 


0.06507 
0.07521 


0.07445 
0.07900 


0.06179 
0.06799 


0.04563 
0.04713 


0.06276 
0.06422 


0.02922 
0.02517 


0.04862 
0.04712 


E = 0.10076 


0.07394 
0.07343 


0.05838 
0.05820 


0.05052 
0.05004 


0.05617 
0.05633 


0.05246 
0.05270 


0.05981 
0.06028 


0.03972 
0.03997 


E = 0.19028 


0.09117 
0.09156 


0.07269 
0.07320 


0.09783 
0.09833 


0.09281 
0.09336 


0.08522 
0.08561 


0.09264 
0.09330 


0.08767 
0.08849 


E = 0.30482 


0.13357 
0.13349 


0.13348 
0.13410 


0.13017 
0.13075 


0.13751 
0.13803 


0.13876 
0.13888 


0.13665 
0.13655 


0.13435 
0.13394 


E = 0.43527 


0.19660 
0.19844 


0.20009 
0.20126 


0.20985 
0.21099 


0.20906 
0.21011 


0.22112 
0.22164 


0.20315 
0.20313 


0.20032 
0.19922 


E = 0.57202 


0.28580 
0.28877 


0.25783 
0.25994 


0.29550 
0.29795 


0.29339 
0.29540 


0.28875 
0.29021 


0.28975 
0.29024 


0.29050 
0.28895 


E = 0.70604 


0.36801 
0.37857 


0.39740 
0.39802 


0.39422 
0.39841 


0.39359 
0.39638 


0.39232 
0.39354 


0.37328 
0.37355 


0.38446 
0.38460 


E = 0.82974 


0.46159 
0.46487 


0.46709 
0.47145 


0.48083 
0.48596 


0.46971 
0.47334 


0.47193 
0.47344 


0.47730 
0.47800 


0.48521 
0.48630 


E = 0.93767 


0.50948 
0.51818 


0.53934 
0.53885 


0.53530 
0.53656 


0.54047 
0.54306 


0.52785 
0.52760 


0.53862 
0.53854 


0.54295 
0.54123 


£■ = 1.02672 


0.56526 
0.56125 


0.56621 
0.56722 


0.58570 
0.58918 


0.58284 
0.58643 


0.57729 
0.57799 


0.58613 
0.58700 


0.57540 
0.57636 


E = 1.22634 


0.63279 
0.64082 


0.65405 
0.65811 


0.65635 
0.65997 


0.65077 
0.65550 


0.64431 
0.64588 


0.65041 
0.65104 


0.64915 
0.65008 



TABLE I: Maximal Lyapunov exponents Ag ~ aXa (upper entries) and Ah ~ aXn (lower entries) as a function of total energy 
E — g^Ea, Higgs self-coupling R and number TV of lattice sites per dimension. 



passing that chaos is not only damped by gauge-field masses generated via spontaneous symmetry breaking, but also 
by those due to quantum fluctuations according to the Coleman- Weinberg mechanism |34|. topological excitations, 
polarization of the heat bath at finite temperature, and external charges Q.) In Ref. [3^ chaotic behavior was 
observed [ssj only beyond the threshold energy Eth — 0.3 (showing that chaos is not universal in the broken phase), 
and for the energy E = 5.07 ^ Eth the MLE was found to be A ~ 0.25 [1^, i.e. an order of magnitude smaller than 
our value A ~ 2.75 in the unbroken phase (which we linearly extrapolate [cf. Sec. IIV B| from the values in Table 1 up 
to E = 5.07). Since constant fields with their few degrees of freedom can exhibit stronger chaoticity and thus produce 
larger MLEs than our randomized initial configurations, this comparison gives a quantitative idea of how much the 
chaotic YMH instability is damped by the Higgs mechanism in the broken phase. 

Another issue which can be addressed quantitatively on the basis of the data in Table 1 is the relation between the 
maximal Lyapunov exponents Aq and Ah, which are obtained from the gauge and Higgs field distance measures (|29)) 
and respectively. This relation was subject to some debate, in particular at strong coupling [1^, H^l- After an 
exploratory study in Ref. [20] , Ref. ^27} provided a first lattice estimate for YMH theory. The Ah extracted from the 
growth rate of the Higgs field distance measure was found to become smaller than Ag when the Higgs self-coupling R 
increases. At k = 24 and for TV = 10, in particular. Ah was estimated in ReL f27| to be about 15% smaller than Ag- 
Comparison with the static gauge and Higgs boson damping rate in (resummed) thermal perturbation theory then 
cast doubt on their relation to the same Ag and led to the speculation that the Higgs damping rate may instead be 



related to Ah |27[. These ideas were later questioned in Ref. [26| whose improved calculation found Ag and Ah to 



agree, although only for one trajectory pair at fixed energy and S = 0.2, k — 1, R — 0.2 and TV — 10. 

In our case, the Ag (upper entries) and Ah (lower entries) values in Table 1 agree within errors (at the percent 
level) in all of the covered YMH phase space, with the deviations slightly decreasing for increasing E and R. Our 
results therefore show that the finding of Ref. [1^ was not an accidental outcome of one specific initialization choice 
but that indeed 



Ag {E,R,N)~~Xh {E,R,N) 



(41) 
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Since the individual relaxation times r of the gauge and Higgs sectors are set by the inverse MLEs, i.e. tg.h — ^q^h, 
Eq (|¥T|) naturally explains the observation tg ~ th in Sec. IIII Bl i.e. the fact that both gauge and Higgs sectors 
(separately) self-thermalize over about the same relaxation time. Equation (|4ip further squares with the general 
expectation that the maximally unstable field mode of a dynamical system, i.e. the mode associated with the MLE, 
dominates the exponential distance growth. Hence the MLEs should be independent of the metric used to extract 
them (modulo constant factors which depend on the field powers involved in the definition of the metric). A possible 
exception to this rule may arise, however, if the distance measure is blind to the maximally unstable eigenmode. In 
Ref. it was argued that such a situation occurs in YMH theory at large coupling R, where the quartic Higgs self- 
interaction dominates the potential Higgs energy (jl2p : the amplitude i?^. then remains practically unchanged during 
time evolution and decouples from the maximally unstable gauge-field mode to which the Higgs distance measure (I30p 
consequently becomes insensitive. At the still relatively weak coupling R, = 24, and the energy E ^ 0.8 of Ref. [23| 
where dn (t) is rather strongly time dependent [cf. Fig. [5|ay|, however, our results for Ah and Ac differ by only about 
1%. This suggests that the 15% deviation found in Ref. (231 should mainly be attributed to numerical uncertainties. 



B. Energy dependence and relation to plasmon damping rate 

As already mentioned, the dependence of the MLEs on the gauge-field energy per degree of freedom, Eq = Ejj2 (cf. 
Eq. (IM])). and on the total energy E of the YMH system is a particularly important issue. In pure SU(2) 0,11)] and 
SUf3) j2d [ssj Yang-Mills theory (whose scaling properties imply that the dimensionless Lyapunov exponent A = Aa 
can only depend on Ep = g^aEp), the approximately linear relation 

Ajv. ~ CN^g'^Ep (42) 

with C2 — 0.17 and C3 ~ 0.10 was established numerically in the weak-coupling regime. (An improved SU(2) analysis 
and a careful discussion of the involved errors [25^ , triggered by questions raised in Ref. [2y] , later confirmed the results 
of Refs. J,, 2^ ) The empirical relation helps to clarify the physical role of the MLEs in hot quantum gauge theory. 
Since the Lyapunov exponents were extracted at times < tq ~ A~^, i.e. after the gauge sector has preequilibrated 
[83 |. the thermal gauge- field ensemble has according to Eq. (p7l) reached the temperature Tq = E.p/2 (at sufficiently 
large average plaquette energy E-p). Together with Eq. (j42p this implies the linear relationship Xn^ — ^n^Tq, and 
comparison with the static plasmon damping rate 70, Af^ of hot quantum SU(A'^c) YM theory, as calculated to leading 
order in hard-thermal- loop rcsummed perturbation theory [5|, then revealed the at first rather unexpected relation 

Ajv. ~ (43) 

for Nc — 2, 3. [The factor of 2 arises because the growth rate of the distance ((29|) is twice that of the distance between 
the gauge fields.] Subsequently, Eq. (j43p has been derived under a few heuristic assumptions (in particular on the 
ergodicity of the gauge-field evolution) in Ref. ,9] . 

On the basis of the rather exhaustive data set in Table 1, we are now able to address the analogous question of how 
the MLEs are related to the average plaquette and total energies in the weak-coupling regime of YMH theory. In Fig. 
ini we plot the MLEs for k = 1 on lattices with N G {6, 10, 20, 30} (corresponding to the first four columns of Table 1) 
in the full range of average plaquette energies < ag'^Ep = Ep < A. Figure [7] contains all remaining MLEs of Table 
1, i.e. those for R = {1, 8, 16, 24} at = 20. The straight lines also drawn in Figs. [S]and[7]are the best linear fits to 
the data: 

Ag,h (Ep) ~ O.Uag^Ep ~ ^Ep. (44) 

The figures show that the MLEs indeed depend within errors linearly on the average energy Ep per plaquette, as in 
YM theory. In fact, the linearity of A (Ep) seems to be a nontrivial consequence of the non-Abelian nature of the 
gauge group. (The MLEs of scalar theory and Abelian U(l) gauge theory, in contrast, were found to vanish in 
the continuum limit a ^ I20| .) Remarkably, even the slope of the linear relation (|44l) is almost identical to that in 
SU(2) Yang-Mills theory [tI I20II2H. (It is also consistent with the value of the ratio X/Ep which was extracted from 
the trajectory with S = 0.2, k = 1, ^ = 0.2 and A^ = 10 in Ref. 

Equation implies that for identical gauge field energy the MLEs of YM and YMH theory are approximately 
equal. This provides our main evidence for the maximally unstable YMH mode to belong primarily to the gauge 
sector, and suggests that the chaoticity and equilibration properties of the Higgs sector are mediated by this gauge- 
field mode as well (at least at weak coupling and if the major part of the initial energy is stored in the gauge sector). 
It also makes it more plausible that the MLEs of YMH theory are related to the gauge field damping rates [1^ |2g] . 
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Furthermore, it is consistent with the approximately equal relaxation times tq th (cf. Sec. IIII Bl) and exponential 
distance growth rates, cf. Eq. (|1T|) . in the gauge and Higgs sectors. 

Nonetheless, the plaquette energy dependence of the MLEs in Figs. [6] and [7] also shows small systematic deviations 
from linearity which become most notable towards the lowest Ep values. The same effect was observed in pure Yang- 
Mills theory ^] , and a glance at the criterion ((39|) indicates that finite-size errors are responsible for the systematic 
upward trend of the MLEs at the smallest Ep. In fact, this is what one would intuitively expect since field modes 
with longer average wavelengths are more strongly deformed by the periodic boundary conditions. The slopes of the 
logarithmic distance histories become most strongly modulated towards smaller Ep (cf. Sec. IIII C[) . furthermore, 
which introduces additional systematic errors. Together with the finite-time errors to be discussed in Sec. llVCI thev 
might cause additional deviations from a linear energy dependence of the MLEs. Towards the maximal value i?p = 4 of 
the average energy per plaquette, on the other hand, lattice spacing [cf. Eq. ([55)1 ] and compact phase-space artefacts 
are likely to affect the results [l^, [s^ ■ 

Since the YMH system has a second characteristic energy scale besides Ep, i.e. the total energy E which additionally 
includes both the energy stored in the Higgs field and in the gauge-Higgs interactions (cf. Sec. IIII Bp and is strictly 
conserved at all times, it is natural to ask how the MLEs depend on E. In order to answer this question, we plot our 
MLEs in Figs. [5] and [5] as a function of E (for the same R and N values as in Figs. [S]and[7]) and find the dependence 
on the total YMH energy to be approximately linear as well: 

Ag,h {E) ~ 0.55E. (45) 

The above scaling behavior can be understood by recalling that our MLEs were extracted during evolution times 
t/a < 400 over which the distances generally saturate, but before the gauge and Higgs fields have started to exchange 
appreciable amounts of energy. A glance at Fig. [2] shows that after the gauge and Higgs sector have separately 
preequilibrated (i.e. for t/a 3> A~^), Eq = Ep/2 and En are practically time-independent in this phase. Moreover, as 
mentioned in Sec. IIII CI staying in the weak-coupling regime requires initial conditions which (except for the smallest 

5 = 0.05 ) store considerably more energy in the gauge than in the Higgs sector (cf. e.g. Fig. [2] which corresponds to 

6 — 0.2). In this situation one derives from the definition of E in Sec. IIIIBI which implies Ep = —AEn/S + lOE/3, 
and from Eq. (gi]) that 

Ag,h {E,En) ^ -^^H + Ie ^E, (46) 

which explains the linear behavior and numerical slope of Eq. [Eq. also explains the numerical scaling 

relation Ac (e) ~ 0.53E' 24] for SU(2) YM theory where Eh = 0.] We reemphasize that these results hold for MLEs 
extracted in the time window A~^ -^t/a < 400 during which Eq and i?H <C -Eg remain practically constant. 

C. Long-time evolution of the Lyapunov histories 

At later times the Higgs sector will pick up energy from the gauge sector, i.e. Ep will drop (for 5 > 0.05 ) while E 
remains constant (cf. Figs. [2] and [3]). In the t ^ oo limit the MLEs must attain a constant value, as implied in their 
formal definition, and so will Ep. This saturation is strongly delayed, however, by the exceptionally long relaxation 
times which govern the equilibration between the gauge and Higgs fields. In the remainder of this section we will 
analyze the quantitative impact of this saturation behavior on the extracted MLE values. To this end, we compute 
the "Lyapunov histories" 

Ag3 (E, -.,N;t):=^ In '^^^'^l = ^ E In Ag,h {E, N) (47) 

t dl j [E, K, N; 0) t 

[where (Igh^ are the rescaled distances (|3T]) and ((32)) . and the Sk are the rescaling factors obtained after the fc-th 
scaling step with rescaling period r], which approach the exact MLEs in the t ^ oo limit, for eight long-time field-pair 
trajectories in the time interval t/a E [0,20000] on an = 10 lattice. (In order to improve numerical efficiency, 
we increase the rescaling period r with increasing saturation time, i.e. with decreasing d, as detailed in the figure 
captions below.) 

In Fig. [To] we show the Lyapunov histories Ag (i) /E (black lines) and Ah {t) /E (grey Hues), "normalized" by 
the total energy, for the four long-time field-pair histories with R = 1 and initial magnetic energies specified by 
6 = 0.2,0.3,0.45 and 0.6. (This corresponds to approximately equally spaced Ep values, cf. Figure [1]) A first 
important characteristic of all Lyapunov histories is their monotonic decrease with time. Moreover, the saturation 
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of Ag (t) for large t can be seen to proceed very slowly: especially for smaller 5 it is not fully completed even at 
t = 2x lO^a. In all four cases, furthermore, Aq (t) starts out somewhat larger than Ah (t) but becomes smaller when 
the gauge and Higgs sectors start to exchange substantial amounts of energy. The deviations between Ag (t) and 
Ah (t) remain at the one-percent level during the initial time evolution (as reflected in the Xq and Ah estimates of 
Table 1) and increase systematically up to 5% at t = 20000a. Hence Ag (t) and Ah {t) remain approximately equal 
over the whole time evolution until they turn into the MLEs for t —^ oo. 

In Fig. [Tl]we display the Lyapunov histories Ag,h (t) /E for an intermediate S = 0.3 and the four values R = 8, 16, 24 
of the Higgs self-coupling. The main tendencies observed in Fig. [TO] remain intact for larger k, although increasing 
Higgs couplings further delay the saturation of the Lyapunov histories. Indeed, already for k = 8 it is more difficult 
to reliably extrapolate Ag,h {t) from the simulation interval t/a G [0, 20000] to the MLE value at t ^ oo. On the 
other hand, larger k values further reduce the deviations between Ag (t) and Ah (t) (which suggests that the opposite 
tendency observed in Ref. [27| was due to a numerical artefact), and they also reduce the long-time variations of 
Ag,h (t) /E and hence the finite-time errors of the MLEs. 

In order to check whether the used integration time step At = 0.001 is small enough, we have also performed a 
simulation with half of its value for R, ~ 24. The corresponding curve, also drawn in Fig. 111! is essentially identical to 
the one with the larger time step, which shows that the latter has no relevant time discretization error. Finally, we 
note that when the Lyapunov histories decrease during equilibration, one may expect their sensitivity to the Higgs 
sector to become larger. Figure [TT] shows that their dependence on the Higgs self-coupling R is negligible at smaller 
evolution times (t < 2000a), as manifest in Table 1, but indeed becomes more pronounced for larger t. 

We now turn to the examination of the ratios Ag (t) /Ep {t) and Ah (t) /Ep (t) which we plot in Figs. [T^ and [TOl 
for the same parameter values as in Figs. [TOl and [TT] The closely parallel movement of Ag (t) and Ah (i), and the 
systematics of their small deviations remain visible here as well. The initial drop in Ag,h (t) / Ep (t) (in particular 
for 6 = 0.2) falls into the time period during which Ep (t) is practically constant, i.e. it is caused by the decrease 
of Ag (t) ■ Later on the gauge-field energy Ep (t) starts to drop (cf. Fig. [2]) and overcompensates the continuing 
decrease of Ag,h (t). This causes the ratios Ag,h (t) /Ep (t) to rise. For 6 = 0.2, one furthermore finds from (shghtly 
extrapolating) Fig. [5] that Ep (t) = 2Eq (t) 2E for t > 13000a, so that the continuing, slight decrease of 

Ag,H jt) t/a>13000, 5=0.2 1 - 

-g^ ^^G,Hit) (48) 

for t > 13000a has again to be attributed solely to the behavior of Ag,h (t)- An important result, visible in both 
Figs. [12] and [TO] is that the ratios Ag,h {t) /Ep (t) saturate significantly earlier than Ag,h (t) /E even at larger values 
of R. [Short-time fluctuations of the average plaquette energy Ep (t) cause the time evolution of Ag,h (t) /Ep (t) to 
appear more ragged.] This indicates that for large t, i.e. on the approach to full equilibrium, the average gauge energy 
Ep (t) /2 decreases at the same rate as Ag,h {t). 

As already alluded to, our long-time evolution results allow for a quantitative assessment of the finite-time errors 
in the MLE estimates of Table 1, which were extracted at rather short evolution times. Figure [TO] indicates that for 
R = 1 the time variations of Ag (t) reach about 25% for 6 = 0.2 and about 30% for 6 = 0.6, while they remain about 
5% smaller for the corresponding Ah (t). These variations may be considered as a (conservative) upper bound on the 
systematic finite-time errors, in particular for larger S and smaller R values, and on the corresponding overestimates 
of the Ag,h in Table 1. (The tendency to overestimate the MLEs when extracting them at shorter evolution times 
was also noted in Refs. [Slli.) 

Finally, our long-time analysis allows us to clarify what happens to the two scaling laws ([TT]) and ([T5]) on the approach 
to total equilibrium in the t oo limit where the Lyapunov histories saturate. In fact, a rather reliable extrapolation 
to this limit (in particular for larger energies) can be achieved by taking advantage of empirical evidence for the 
asymptotic evolution-time dependence ~ with which the Lyapunov histories approach the MLEs [2^. (This 

finite evolution-time behavior is analogous to the finite-size behavior ~ {Na)~^^'^ [H, which results from sampling 
ergodic states (23j . On similar grounds, the width of Gaussian fiuctuations around the mean value of the average 
Lyapunov exponent was argued to decay as HI]-) Hence the new functions Ag,h can approximately 

be fitted by straight lines (with a potential systematic bias) and the MLEs determined as the intersections with the 
t~^/^ = axis [25I]. Extrapolating the curves in Figs. [TO] and [TT] in this way to infinite evolution time yields our best 
estimates for the MLEs which we plot in Fig. [14] as a function of Ep (upper panel) and E (lower panel). These figures 
show that both dependencies remain to good accuracy linear. The best linear fits (also shown in the figures) are 

Ag (^p) =^ 0.149^p, Ah (^p) ~ 0.163^p, (49) 

and 

Ag {e) ~ 0.308i?, Ah {e) ~ 0.338£;. (50) 
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Equations ([49]) show that to an accuracy of at least about 10% the scahng law (|44)) found before full equilibration, 
with the same coefficient as in pure YM theory, indeed remains intact asymptotically. (First indications for this 
behavior were observed in Ref. [2y|_on the basis of a trajectory for t < 12000a.) Qualitatively, this is also reflected in 
Fig. [T^ where the ratios Ag,h {t) /Ep (t) start at around 1/6 and asymptotically return to it for very large t (while in 
the meantime deviating by maximally (i.e. for the largest 6) about 20%, mainly when most of the energy is exchanged 
between the gauge and Higgs fields). 

The Eqs. (HH) also explain the E dependence of the asymptotic Lyapunov histories exhibited by the fits (|50p . 
Indeed, after full equilibration at t ^ oo with Eq = Ep/2 = En = E one expects from Aq.h (^^p) = acH-^p — Ep/6 
that 

Ag,h {E) = 2aG,HE ~ ^E (51) 

which is within errors identical to Eqs. ([50|) . Hence in equilibrium the linear dependence of the MLEs on Ep implies 
a linear dependence on E with twice the slope (as realized to good accuracy in the fits (|49l) and This fact 

went probably unnoticed in Ref. [1^ which argued against linear scaling of the Ag,h with E on the basis of field 
evolution trajectories over maximally several thousand lattice time units, i.e. likely too short to bring the system 
close enough to equilibrium. (In addition, our above observation supports the diagnosis of Ref. [13] which attributes 
the logarithmic energy dependence found numerically for several SU(2)-YM MLEs after very long evolution times 
to finite-size artefacts of the monodromy-matrix method.) More generally, linear scaling of Ag,h with Ep implies a 
linear dependence on E in any time window during which Ep cx E. This condition seems to be satisfied only when 
the gauge and Higgs sectors do not exchange relevant amounts of energy, however, i.e. only in the preequilibration 
phase and after mutual equilibration is essentially achieved. Nevertheless, the slopes of the E dependence in these 
two time intervals are different (5/9 and 1/3, respectively). 

The above evidence for the linear dependence of the MLEs on the average magnetic energy Ep to prevail for 
t —^ oo appears consistent with our previous indications for the maximally chaotic mode to reside mainly in the gauge 
sector, with our finding that the scaling behavior (|^^ sets in way before the gauge fields have full access to the energy 
stored in the Higgs sector, and with the result that the ratios Ag,h (t) / Ep it) saturate significantly faster than the 
Ag,h it) themselves. 

D. Mciximal Lyapunov exponents of initially homogeneous magnetic fields 

In this section we digress from our main subject and apply some of the numerical techniques developed above to 
the time evolution of spacially constant, non-Abelian magnetic fields. In the pioneering days of QCD such homoge- 
neous magnetic fields were perturbatively established to be unstable in pure YM theory [37] . This so-called Savvidy 
instability was later explored in the nonperturbative domain by numerical methods similar to ours [2, i38| . It provided 
early indications for the complexity of the Yang-Mills vacuum and has triggered the development of stochastic and 
chaotic concepts for vacuum structure and quark confinement [TtI . [39| . In the following we are going to study the 
impact of the matter (i.e. Higgs) fields on the Savvidy instability. 

As a benchmark for comparison with the YMH case, we first reproduce the nonperturbative time evolution of the 
distance between initially adjacent, homogeneous magnetic fields in pure YM theory on an iV = 10 lattice. The 
non-Abelian magnetic field is defined as Bp — arccosir [Up), and the fields are initialized with total energy E = 0.57 
by setting E'^^^ (0) = 0, S^,i2 (0) = 0.899, B^,i3 (0) = 0.791 and B^,23 (0) = 1.453 for aU x. I n Fig. [I5l we compare 
their logarithmic distance evolution to that of an initially randomized gauge field (cf. Sec. IIII Ap with the same 
energy. The constant magnetic field pair turns out to have about twice the average slope of In (Iq (t) in the linear 
region, i.e. the homogeneous magnetic field is substantially less stable than the random field. This result corroborates 
similar findings in Ref. [3]. 

We now turn to the analogous time evolution of initially constant magnetic fields in YMH theory, again on a lattice 
with iV = 10 sites per dimension. As in all previous sections, the Higgs field is initialized at the spacially constant 
value Rx {0) — R = 0.2, and the initial values 4>{0) = 0, E"^ j (0) — are imposed in order to satisfy Gauss' law 
([27)1 . The B field is initiafized at the values ^^42 (0) = 2.319,' S^, 13 (0) = 2152 and Br, 23 (0) = 1.428 for aU x. We 
further choose S = 0.3 and k = 1 in order to inject the same total energy E = 0.57 as in the YM case above. The 
logarithmic distance evolution under these conditions is displayed in Fig. [12] for the gauge and Higgs field metrics P5)) 
and again together with its counterpart for a corresponding random field. As in the YM case, the homogeneous 
magnetic field produces about twice the slope in the linear region. Hence the presence of the matter fields seems 
neither to dampen nor to enhance the increased instability of the homogeneous magnetic field configurations relative 
to a random configuration of the same total energy. This seems to be consistent with our above evidence for the 
chaoticity of YMH theory to be dominated by the gauge sector. 
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The slopes of In da (t) for both constant magnetic and random fields, however, are (in the linear region) about 
twice as large in the YMH example of Fig. [11] than in YM theory (Fig. [TS]) . Probably this result depends rather 
strongly on the initial conditions, and especially on how they distribute the initial energy over the gauge and Higgs 
field sectors. Our above initial conditions were chosen to provide a demonstrative example for the presence of matter 
fields to strongly enhance the Savvidy instability. This massive impact raises the possibility that simulation results for 
those gauge-field instabilities which drive isotropization and thermalization in the aftermath of high-energy nuclear 
collisions could be significantly modified in the presence of quark fields as well (although they have smaller occupation 
numbers). 

V. RELATED EQUILIBRATION PROCESSES IN COSMOLOGY AND NUCLEAR COLLISIONS 

In the following section we are going to discuss several aspects of nonequilibrium processes in the early Universe and 
in the aftermath of high-energy nuclear collisions which are pertinent in our context. We comment on the impact of the 
chaotic thermalization properties calculated above and on results of classical gauge-theory simulations related to ours. 
We also suggest a few promising extensions of our work which would help to clarify the role of chaotic thermalization 
processes after nuclear collisions and the nature of the "apparent" or "pre-" equilibrium at the beginning of the 
subsequent hydrodynamic evolution phase. 

A. Early Universe 

According to the inflation paradigm, the vacuum energy of one or more classical, scalar inflaton fields dominated the 
very early Universe. This energy caused a typically exponentially accelerated expansion period which very efficiently 
diluted particles and fluctuations [ioj . It left the universe in a supercooled, highly nonthermal state which was 
practically devoid of matter, radiation and entropy. Hundreds of models for the phenomenologically very successful 
inflationary scenario were proposed A compelling "microscopic" theory, however, which is able to explain e.g. 
the nature of the inflaton(s) and the very specific properties of their dynamics, has not yet been established, partly 
because it involves unknown physics beyond the standard model. 

For the post-inflationary reheating period, during which the Universe thermalized at a still very large "reheating 
temperature" , the theoretical situation is similar: there exist many scenarios and model calculations for specific 
processes whereas the underlying dynamics as a whole is not yet settled. During reheating a huge amount of entropy 
was released. All the matter and radiation of the present Universe was created and the energy density of the 
inflaton(s) was transformed into a hot and ultrarelativistic plasma. (Afterwards the Universe expanded essentially in 
the Friedman- Robertson- Walker geometry and cooled almost isoentropically according to the "hot-big-bang" scenario.) 
In the following we will be particularly interested in the (semi-) classical phases during the reheating period, when the 
large occupation numbers of the participating field modes made contributions from chaotic thermalization relevant. 

The probably most important of these phases, referred to as preheating, is suggested to have taken place immediately 
after inflation [4l| . During this very short period, which lasted about 10~^^ sees., particle production became explosive. 
Preheating can be induced either by a tachyonic instability of the inhomogeneous modes which accompany electroweak 
symmetry breaking |42|], or by the stimulated decay of an almost homogeneous inflaton which coherently oscillates 
with an initial amplitude of the order of the Planck mass. In the latter case, the accelerated decay is the consequence 
of a parametric resonance with condensates composed of the produced particles [4T| . A detailed understanding of the 
preheating process is particularly crucial because the bulk of the initial conditions for the subsequent thermal history 
of the Universe are settled at its end. Since the homogeneous energy density of the inflaton transfers exponentially 
rapidly into highly occupied, inhomogeneous out-of-equilibrium modes, furthermore, the violently nonperturbative 
and nonequilibrium processes underlying preheating are amenable to classical lattice simulations [43] . 

After preheating, a variety of crucial thermalization processes began to drive the Universe towards equilibrium. 
Classical lattice simulations of such processes indicate, furthermore, that the infrared modes excited during preheating 
evolve towards a saturated occupation number distribution long before thermalization completes |44|. Such effects 
have been interpreted as signs of "prethermalization" , characterized by an energy-pressure relation approximating 
an equation of state [H, . Later on, the distributions move towards complete saturation by cascading towards 
ultraviolet and infrared modes (as in Kolmogorov wave turbulence). During the last and longest stage of equilibration, 
finally, the particle distributions become fully thermal ^8^. Simultaneously, the occupation numbers drop until 
quantum physics eventually dominates and classical simulations become ineffective. 

Among the many important nonequilibrium processes which shape the reheating period are backreactions on the 
inflaton field which eventually stop particle production, violent and still nonperturbative particle rescattering events 
which very efficiently generate entropy, the nonthermal production of heavy particles as well as phase transitions. 
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Further reheating processes which left crucial imprints in our present Universe are primordial magnetic field generation 
p^ . topological and large-scale structure formation as well as baryogenesis 11 1, which created the observed abundance 
of baryons over antibaryons. In our context, baryogenesis is particularly interesting since it can efficiently proceed 
only far from equilibrium. As already mentioned, lattice simulations of classical real-time field evolution are a method 
of choice for the analysis of such processes. In fact, they turned out to be particularly useful for shedding light on the 
preheating phase whose immense particle production rate and condensate formation require a fully nonperturbative 
treatment while the generated, large (boson) occupation numbers ensure a classical field evolution [47l |. 

In the following, we will focus on classical lattice simulations of cosmological pre- and reheating processes which 
were based on dynamics including the gauge-scalar sector of the standard model [86[, i.e. the SU(2) YMH theory 
with a fundamental Higgs doublet ([T]) [83|- This theory has been used, notably, to explore the electroweak symmetry 
breaking transition [48|. It underlies our own work, furthermore, and thus allows for some partial and qualitative 
comparisons - although the initial conditions (and the additional field content including e.g. an inflaton) required to 
describe cosmological situations differ considerably from the random ones which we have adopted above. 

The analysis of electroweak baryogenesis at energies of the order of 100 GeV is, as already alluded to, an especially 
interesting application of the YMH model ^ and its extensions jjl!, ^49*1 . Such studies are mainly motivated by the 
question whether "minimal extensions" of the electroweak standard model, which implement e.g. additional neutral 
scalar inflaton field(s) or CP violating couplings, are able to explain the observed baryon asymmetry of the present 
Universe. As a step towards clarifying this issue, the baryon (and lepton) number and CP violating sphalcron tran- 
sition rate and the Chern-Simons number diffusion in the unbroken phase were studied in Refs. [15, 50, 51, 52]. In 
a hybrid inflation scenario [s^, based on an additional singlet inflaton which couples to the Higgs field to accelerate 
particle production (compared to only gravitational coupling), the nonequilibrium preheating dynamics was found 
to generate Chern-Simons number, a prerequisite for electroweak baryogenesis, locally and stochastically [s^ ]. (For 
some cautionary remarks on the reliability of classical lattice results in this preheating scenario for baryogenesis see 
Ref. [5^.) In the same dynamical framework, primordial magnetic fields are produced with sufficient magnitude and 
correlations to act as seeds for the magnetic fields observed in galaxies and galaxy clusters today [i^. Electroweak 
baryogenesis during a cold electroweak transition with tachyonic preheating (induced by a spinodal Higgs field insta- 
bility) and additional CP violation generated by a coupling of the Higgs field to the topological charge density of the 
gauge field, has been investigated in Ref. [5^. The particle distribution functions of the Higgs and gauge fields (on 
which kinetic theory is based) were extracted in Ref. [575 from the correlators of the simulated classical fields, with 
the electroweak phase transition modeled by a quench. 

Our main goal in the present paper was to study generic chaotic thermalization properties of YMH theory on the 
basis of random initial conditions. In this respect, our work differs from the more specialized simulations discussed 
above (which partly also include additional dynamics). Although this prevents a quantitative comparison of the 
results, we believe that the chaotic thermalization processes which we have analyzed should be relevant for most of 
the mentioned post-inflationary pre- and reheating processes as well. By adapting the initial conditions to cosmological 
situations, in particular, one could directly investigate the contributions of deterministic chaos, as measured by the 
Lyapunov exponents, to specific nonequilibrium processes. It would for example be interesting to follow the sphaleron 
transition and magnetic field production rates during the different stages of chaotic thermalization. It would also 
be useful to study the impact of the lattice spacing on the chaotic thermalization rates [ssj and to include physics 
beyond the standard model, e.g. an inflaton fleld, into the analysis. 



B. Nuclear collisions 



Thermalization properties of excited quark-gluon matter, as produced at the SPS, RHIC, LHC and (future) FAIR 
colliders, have been intensely studied in various approaches of increasing sophistication [88| (see Refs. [6^ for recent 
reviews). The detailed local equilibration mechanisms are of central importance for the heavy- ion programs since 
they determine whether a new, deconfined state of matter - the quark-gluon plasma - can be locally thermalized 
in the aftermath of high-energy nuclear collisions, i.e. whether the produced system equilibrates fast enough for 
thermodynamic concepts to apply before it disintegrates. 

The thermalization issue became even more intriguing when RHIC results showed that the produced matter starts 
to behave collectively after times of less than 1 fm/c and is subsequently described by essentially ideal Bjorken 
hydrodynamics with almost maximal elliptic flow 61]. These findings are generally interpreted as a surprisingly fast 
apparent thermalization of the system, which minimally requires the isotropization of the long-wavelength modes 
participating in the hydrodynamic behavior 62] and perhaps the onset of prethermalization [45]. In any case, the 
very short (pre-) equilibration time cannot be explained by weakly coupled parton-parton collisions alone (60l. [63j. 

The time evolution of a typical RHIC reaction (which is sometimes referred to as a "little bang" to emphasize 
similarities with the big bang of the Universe) begins with very hard initial interactions between the high-momentum 
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partons of the colliding nuclei. These generate the highest-momentum particles in the final state. Afterwards, at 
about t ~ 0.2 fm/c, most of the soft particles in the final state are coherently produced and form a nonequilibrium 
system (sometimes called a "glasma" [i^]) of high energy density. The "bottom-up" thermalization scenario (65l | 
assumes that the initial properties of this system are determined by the QCD saturation mechanism (soj . i.e. that 
they are dominated by coherent small-x gluons of a very high density. These gluons originate from the low- a:; part 
of the nuclear wavefunctions (x = pt/ \/s where s is the total energy) and carry transverse momenta pt of the order 
of the saturation scale Qs ■ Since for RHIC collisions Qs ^ i GeV ^ Aqcd , this initial state should be amenable to 
weak-coupling techniques [lOl- In the color glass condensate (CGC) model 68], for example, the highly populated 
small- a; gluon states are treated as the soft modes of classical Yang-Mills fields with typically large amplitudes while 
the hard field modes are represented by static sources. 

More generally, as long as the gluon mode occupation numbers stay large enough to suppress quantum effects, 
their evolution can be described in terms of classical gauge fields which may be simulated fully nonperturbatively 
on a spacial lattice. Over the last years, an increasing amount of such numerical simulations was performed for the 
gauge groups SU(2) and SU(3) in one, two and three spacial dimensions. The dynamical settings included the CGC 
model [69|, , hot-thermal- loop (HTL) effective theories (equivalent to a coUisionless Vlasov equation) [H, [tiI, [T^I 
and equations of Wong- Yang-Mills type [tI, I?!, [tHIj subject to minijet or bottom- up initial conditions. Among the 
calculated observables were e.g. energy densities and gluon multiplicity distributions. 

Numerical simulations confirmed, in particular, that sufficiently anisotropic parton momentum distributions, as 
typically produced in heavy-ion collisions, can induce the onset of a very fast, collective isotropization and (pre-) 
thermalization [9i| driven by non-Abelian filamentation instabilities [73| . These Weibel-type plasma instabilities 

Pierate an initially exponential growth of the soft-mode occupation numbers even at relatively weak couplings 
, H^, [zE|- The expansion of the system can strongly reduce this growth, however, and the non-Abelian self- 
interactions appear to rapidly limit it (in three spacial dimensions) to at most linear growth [73 |. This happens 
when the (magnetic) energy deposited in the soft fields returns sufficiently fast to the hard fields, either by plasmon 
excitation effects similar to Kolmogorov wave turbulence [tI, [T^I or via a rapid avalanche [T^I ■ Applications of the 
CGC and HTL effective theories are limited to the weak-coupling and small-amplitude regimes, however, which ensure 
a sufficient scale separation between the hard parton and soft gauge- field momenta [92i] . This technical limitation can 
be avoided by describing the whole system as a classical statistical Yang-Mills ensemble with an UV lattice cutoff 
to substitute for the quantum mechanical suppression of thermal short-distance effects. Hard and soft modes are 
thus originating from the same gauge field and treated on the same footing, which allows for nonperturbatively large 
amplitudes. 

The numerical simulations described in the present paper are based on the same interpretation of the classical 
lattice fields. Although we focused on the thermalization of gauge-Higgs matter by chaotic instabilities, which at first 
seems to be a rather different mechanism, there are several commonalities with equilibration via mean-field plasma 
instabilities. In fact, both are collective processes which shape the momentum distribution of the classical gauge fields, 
and both lead to substantially faster gluon equilibration rates than collisional thermalization. This becomes explicit 
in relaxation times r ~ ~ 0.5 fm/c at typical RHIC energies. (The chaotic relaxation time tends to become even 
smaller for SU(3) gauge fields [11]). As pointed out in Ref. [23|, incidentally, relaxation times of this order imply 
that fluctuations around the mean transverse momentum produced in nuclear collisions are very small [93| (of order 
10~^, i.e. at the percent level, for typical reaction volumes (5 fm)^ and times (0.5 fni) at RHIC), as indeed observed 
in event-by-event fluctuations. 

On the other hand, there are remarkable differences between chaotic and plasma instabilities. The maximally 
chaotic modes are very efficient in generating entropy directly, for instance, while the filament instabilities lead to 
a reversible isotropization (at the mean-field level) which just creates more efficient conditions for the subsequent 
entropy production. Hence Weibel-type plasma instabilities per se seem hardly to affect the MLE values, although 
they could indirectly contribute to our chaotic thermalization if the random initial conditions generate a sufficiently 
anisotropic momentum distribution of the gauge fields. It is tempting to speculate, then, that even under RHIC 
initial conditions the maximally chaotic modes may lead to faster thermalization than plasma instabilities because 
they generate entropy in a probably more efficient one-step process. Moreover, the maximal chaotic instabilities select 
at any time the most unstable direction in phase space and thus remain optimally fast during the entire (classical) 
thermalization process, whereas the filament instabilities are damped by the gauge field's non-Abelian self-interactions 

As mentioned above, a long-term perspective of our work is to understand the role of chaotic thermalization 
mechanisms during the classical evolution phases in the aftermath of nuclear collisions. The focus of the present 
paper was more limited and preparatory, however, namely to map out generic chaotic thermalization properties of 
non-Abelian gauge systems in the presence of scalar matter. (Classical Higgs fields were included because of their role 
in the early Universe and because they can be treated on the same footing as the gauge fields. Quarks, on the other 
hand, would have to be implemented either as particles or as solutions of the Dirac equation in the background of the 
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classical fields, as e.g. in soliton models for baryons [78|.) For this reason, we have chosen random initial conditions 
instead of the more specialized ones required to describe the physical situation after a nuclear collision (or during 
the reheating phase of the early Universe) . Our program could naturally be extended by adapting it to those initial 
conditions with a strongly anisotropic momentum distribution which characterize the hot system created by a high- 
energy nuclear collision. The chaotic properties of plasma- instability-, collision- and cascade-driven thermalization 
processes could then be quantified in terms of the corresponding maximal Lyapunov exponents, thereby relating the 
chaotic thermalization and entropy production rates to the time scales of more conventional equilibration processes. 
(The large difference between the Lyapunov exponents of randomly and coherently initialized fields (cf. Sees. IIV Al 
and II V D|) show that they can depend strongly on the initial conditions.) 

In addition, one may extract a more detailed picture of the chaotic equilibration processes by following e.g. the 
evolution of the field modes' momentum distributions and their anisotropy, as well as the evolution of pressure, 
entropy etc. simultaneously with the Lyapunov histories. For specific comparisons with the effective CGC and HTL 
dynamics, it would also be interesting to follow the growth rate of the gauge field's Fourier coefficients. This would 
further clarify how far chaotic evolution is consistent with the bottom-up thermalization scenario, according to which 
the copiously produced soft gluons (which initially carry only a small fraction of the total energy) draw energy from 
the hard gluons and thermalize very efficiently. The decay of the remaining hard gluons then reheats the soft gluonic 
background until it enters the hydrodynamic evolution stage. In fact, in reductions of the Yang-Mills dynamics to 
a few degrees of freedom, chaotic thermalization was found to start (under suitable initial conditions) among the 
softer modes as well [t^. Moreover, the chaotic behavior of the soft modes turns out to be driven by the soft-hard 
mode coupling, which provides an efficient mechanism for energy transfer from high frequency modes to a low- energy 
multiparticle final state. (Related findings for the scattering of two classical SU(2) YM and YMH wave packets |95| 
were reported in Ref. [soj.) 



VI. SUMMARY AND CONCLUSIONS 



We have investigated chaotic instability and thermalization properties of classical gauge and matter fields in the 
unbroken phase of SU(2) Yang-Mills-Higgs theory on a spacial lattice. Since equilibration proceeds mainly through 
the most unstable field modes, we have focused on a quantitative survey of the most chaotic time evolution patterns 
in terms of maximal Lyapunov exponents (which measure the logarithmic distance growth rates between initially 
neighboring gauge and Higgs fields). 

A main goal of our investigation was to explore the impact of the fundamental doublet of scalar matter fields on 
the chaotic behavior of the gauge dynamics. Towards this end, we first confirmed and extended previous evidence 
for the Yang-Mills-Higgs system to equilibrate over two drastically different time scales: individually, the rather 
weakly coupled gauge and Higgs sectors reach a preequilibrium phase after only a few lattice time units, whereas their 
mutual equilibration is substantially delayed by the matter fields and takes far longer than 10'' units to complete. 
Accordingly, we have generated two sets of maximal Lyapunov exponents for initially random fields: a larger one 
extracted from the field separation rates at the preequilibration stage, and a smaller one obtained from long-time 
trajectories extrapolated to infinite evolution time and thus to full equilibrium. 

The first set was designed to cover a representative part of the weakly coupled phase space and contains about 
80 Lyapunov exponents. In view of its non-negligible finite-time errors, this set was mainly used to study general 
characteristics of the exponents including their energy, coupling-parameter and lattice-size dependence. We found 
the signs of the whole set to be positive (in contrast to results in the broken phase), which implies that chaos is 
at least approximately universal in the symmetric phase. (Possible exceptions may include small nonergodic niches 
as previously encountered in the Yang-Mills phase space.) Yang-Mills-Higgs theories with gauge groups containing 
SU(2) as a subgroup are therefore chaotic as well. In the energy and coupling ranges where both finite-size and lattice- 
spacing errors should be under control, our Lyapunov exponents on lattices with between 10^ and 30'^ sites were indeed 
found to be within statistical uncertainties identical. In addition, we found the maximal Lyapunov exponents extracted 
during the preequilibrium phase to be almost independent of the Higgs self-coupling. This indicates that the nonlinear 
interactions in the Higgs sector provide a relatively minor contribution to the chaoticity of the system. 

In order to survey the important and previously unexplored asymptotic regions of the Yang-Mills-Higgs phase 
space, we have additionally followed several field evolution trajectories over the exceptionally long periods required 
to approach total equilibrium. In particular, we have investigated the long-time behavior of eight Lyapunov histories, 
i.e. logarithmic separation rates between initially neighboring fields, over 20000 lattice units. After extrapolation to 
infinite evolution times they provide our best estimates for the maximal Lyapunov exponents. All Lyapunov histories 
turn out to decrease monotonically with time and saturate relatively slowly even after the energy has attained almost 
complete equipartition. At larger values of the Higgs self-coupling the Lyapunov histories vary less strongly during time 
evolution while their saturation is further delayed. Moreover, the Lyapunov histories divided by the average plaquette 
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energy saturate faster than the Lyapunov histories themselves. This foreshadows an early onset of their linear sealing 
relation in the equilibrated system (see below). The long-time behavior of the Lyapunov histories provides reliable 
estimates for the finite-time errors of the Lyapunov exponents, furthermore, and quantifies in particular how much 
they are overestimated in the preequilibrium phase. 

The physical interest in Lyapunov exponents of gauge fields originates partly from evidence for their linear relation 
to the plasmon damping rates at weak coupling. This relation relies on the empirically identified linear dependence 
of the maximal Lyapunov exponents on the average magnetic gauge-field energy. We have therefore systematically 
scrutinized the accuracy and validity range of this scaling relation on the basis of our full data set. Both during the 
rather long preequilibrium period, i.e. before gauge and Higgs sectors have exchanged substantial amounts of energy, 
and after full thermalization we have indeed found the Lyapunov histories to scale within errors linearly with the 
average gauge-field energy. More specifically, both during preequilibrium and after complete equilibration our results 
establish the relation Aq.h — £'p/6 between the maximal Lyapunov exponents and the average plaquette energy. 
This relation was previously encountered in pure Yang-Mills theory and seems to be a rather exclusive property 
of non-Abelian gauge theories. It implies that the Lyapunov exponents for any given gauge field energy, extracted 
either from the gauge or Higgs field separation rates, remain within errors independent of the presence of the Higgs 
field. This lends additional credence to the suggested equality between the Lyapunov exponents and twice the static 
plasmon damping rates of quantum Yang-Mills-Higgs theory at high temperatures and weak coupling. In addition, 
we have established that both during the preequilibrium stage and after full thermalization the maximal Lyapunov 
exponents also scale linearly with the total energy. This is in contrast to previous expectations and turns out to be a 
consequence of the proportionality between gauge and total energy during both phases. 

We have furthermore studied how interactions with the Higgs field affect the Savvidy instability of constant non- 
Abelian magnetic fields. As a benchmark, we have first obtained the maximal Lyapunov exponent for an initially 
constant magnetic field in pure Yang-Mills theory and found it about twice as large as that of a randomly initialized 
field under otherwise equal conditions. We have then computed the analogous Lyapunov exponent in Yang-Mills- 
Higgs theory for initially homogeneous magnetic and Higgs fields at the same total energy and found it still to be 
about 2 times larger than for randomized fields. Hence the additional matter fields seem neither to dampen nor to 
enhance the Savvidy instability relative to that of random fields of the same energy. Depending on the initial energy 
distribution between gauge and Higgs fields, however, the presence of matter fields can have a strong impact on the 
absolute magnitude of the magnetic field's instability. In order to demonstrate this, we have provided an example in 
which the matter fields approximately double the maximal Lyapunov exponent of the gauge field. 

In all simulations described above, we found the Lyapunov histories to be within errors of at most a few percent 
independent of the underlying distance measure and of whether they were obtained in the gauge or Higgs field's phase 
space. This independence turns out to hold both during the preequilibration phase and towards full equilibrium, and 
both for initially homogeneous and random fields. Hence it confirms the general expectation that the divergence rates 
of the most unstable modes should be equally measurable by any reasonable metric in field space. The behavior of the 
Lyapunov histories depends somewhat on the stage of the thermalization process, however. Before gauge and Higgs 
fields have exchanged substantial amounts of energy, the divergence rate in the gauge sector turns out to be slightly 
larger than in the Higgs sector, while it becomes marginally smaller during later phases of equilibration. A previous 
estimate of an about 15% smaller logarithmic Higgs field separation rate, based on a single trajectory, was therefore 
probably contaminated by numerical uncertainties. 

Our above results strengthen the evidence for the gauge dynamics to provide the main source of chaotic instability 
in the Yang-Mills-Higgs system. The matter fields, in contrast, seem to play a subordinate role (similar to the 
quark fields shortly after a nuclear collision). The evidence includes the fast preequilibration of the gauge sector, 
the observation that the Higgs sector seems to have little impact on the ratio between the constant magnetic and 
random field Lyapunov exponents, the linear dependence of the maximal Lyapunov exponents on the average gauge- 
field energy alone, the finding that this scaling behavior sets in way before the gauge fields have full access to the 
energy stored in the Higgs sector, and in particular the fact that the values of the maximal Lyapunov exponents for 
a given plaquette energy turn out to be within errors identical to those in pure Yang-Mills theory. Moreover, the 
nonlinear Higgs dynamics seems not to contribute substantially to the maximally unstable mode, at least not in the 
range of relatively weak couplings where lattice artefacts are under control. Indeed, at the preequilibration stage the 
Lyapunov histories turn out to be practically independent of the Higgs self-coupling, and even afterwards the coupling 
dependence remains moderate. 

Nevertheless, we found the scalar matter fields to have a major impact on the thermalization of the gauge system. 
As we have shown, their presence can strongly enhance the absolute divergence rate between neighboring gauge 
fields, both homogeneous and random. Moreover, the maximally chaotic mode has an almost immediate effect on 
the Higgs sector, as witnessed by the fact that its separation rate can be monitored equally well by following the 
distance evolution between neighboring Higgs fields. Finally, the presence of the Higgs fields massively prolongs the 
equilibration of the system as a whole, at least at weak coupling. This qualitative effect may be robust enough to 
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prevail in the case of fermionic matter, and hence be relevant for understanding the equilibration properties of the 
highly excited quark- gluon matter produced in ultrarelativistic nuclear collisions. 
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FIG. 1: The average energy per plaquette Sp as a function of the initiahzation parameter 6. 
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FIG. 2: Time evolution of the energy (per degree of freedom) stored in the gauge field Eg (uppermost initially horizontal line) 
and in the Higgs field Eu (lowermost initially horizontal line) and their sum, i.e. the (conserved) total energy E (horizontal 
line). The initially oscillating line starting at zero is the electric gauge-field energy E^i, the one starting at the total initial 
gauge-field energy is the magnetic (potential) gauge field energy Ema-g- (The initial conditions for the underlying trajectory 
were 5 = 0.2, R = 0.2, k = 1, AT = 10 and At = lO"".) 




FIG. 3: Same as in Fig. (2] but for the trajectory subject to the initial conditions S — 1, R = 0.2, k = 1 and A'' — 10. 
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FIG. 4: The logarithmic distance evolution in the gauge (black) and Higgs (grey) sectors at fixed Higgs self-coupling R — 1 
is plotted for four lattice volumina corresponding to A'^ = 6, 10, 20 and 30. The initial magnetic energy is parametrized by 6. 
Panel (a) corresponds to S = 1.0, 0.5, 0.45, 0.4, 0.35, panel (b) to 5 = 0.3, 0.25, 0.2, 0.15 and panel (c) to the two smallest values 
5 = 0.1,0.05. 



27 




10 20 30 40 50 

t/a 

O i i 




50 100 150 200 

t/a 




100 200 300 400 

t/a 



FIG. 5: The logarithmic distance evolution in the gauge (black) and Higgs (grey) sectors on a A'' = 20 lattice for four different 
Higgs self-couplings k = 1, 8, 16 and 24. The distance trajectories are grouped as in Fig. |4] according to their initial average 
magnetic energy (parametrized by 5): panel (a) contains the curves corresponding to 5 = 1.0,0.5,0.45,0.4,0.35, panel (b) 
corresponds to 5 = 0.3, 0.25, 0.2, 0.15 and panel (c) to the two smallest values 5 — 0.1, 0.05. 
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FIG. 6: Values of the maximal Lyapunov exponents Ag (black symbols) and Ah (grey symbols) as a function of the average 
energy per plaquette Ep for k = 1 and N = 6, 10, 20, 30. 
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FIG. 7: Values of the maximal Lyapunov exponents Ag (black symbols) and Ah (grey symbols) as a function of the average 
energy per plaquette Ep for N = 20 and R = 1,8, 16, 24. 
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FIG. 8: Values of the maximal Lyapunov exponents Ag (black symbols) and Ah (grey symbols) as a function of the total energy 
B for K = 1 and iV = 6, 10, 20, 30. 
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FIG. 9: Values of the maximal Lyapunov exponents Ag (black symbols) and Ah (grey symbols) as a function of the total energy 
EfoT N = 20 and k = 1, 8, 16, 24. 
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FIG. 10: Long-time evolution of the Lyapunov histories Ag (t) /E (black lines) and Ah {t) /E (grey lines) for (i5,t) = (0.6, 25) 
(2nd largest initial values), (0.45,30) (largest initial values), (0.3,55) (2nd lowest initial values) and (0.2,110) (lowest initial 
values) (with At = 0.0005, k = 1, E = 0.2 and N = 10). 
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FIG. 11: Long-time evolution of the Lyapunov histories Ag (t) / E (black lines) and Ah {t) /E (grey lines) for k = 8, 16, 24 
(from bottom to top) with At = 0.001. The uppermost curve, overlapping with its larger time-step counterpart, is for k = 24 
with At = 0.0005. (All curves correspond to <5 = 0.3, r = 55, = 0.2 and N = 1Q). 
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FIG. 12: Long-time evolution of the Lyapunov histories Ag {t) /Ep (black lines) and Ah (t) /Ep (grey lines) for 6 = 0.6 (largest 
starting values), 0.45, 0.3 and 0.2 (lowest starting values) (with it = 1 and N = 10). 
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FIG. 13: Long-time evolution of the Lyapunov histories Ag {t) /Ep (black lines) and Ah {t) /Ep (grey lines) for k = 8, 16, 24 
and At = 0.001 (from top to bottom) and k = 24 with At = 0.0005 {S = 0.3, t = 55, R = 0.2 and iV = 10). 
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FIG. 14: Lyapunov histories Ag (t) (black dots) and Ah (t) (grey dots), extrapolated to infinite evolution times, as a function 
of Ep (panel a) and E (panel b). The straight lines are the best linear fits. (For 5 = 0.6, 0.45, 0.3 and 0.2, with R = 1 and 
N = 10.) 
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FIG. 15: Logarithmic distance evolution for initially neighboring, homogeneous, non-Abelian magnetic fields (ragged curve) and 
for initially uoigliboring randomized gauge-field configurations (smooth curve) in YM theory (with N = 1Q and E = 0.572023). 
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FIG. 16: Evolution ol the logarithmic distances (|29[) (black lines) and (|30|l (grey lines) for two initially neighboring, homo- 
geneous, non-Abelian magnetic fields (ragged curves with larger slopes in the linear region) and for two initially neighboring 
random gauge-field configurations in YMH theory (with S = 0.3, R = 0.2, k = 1, = 10, ^ = 0.572023). 



